% m_welfare.m 

clear all; 
close all;
format shortG;
% ==================== LOAD FUNDAMENTAL PARAMETERS ========================
load parameters.mat;                            
load calib_productivity_amenity.mat;     
load calib_setup.mat;   
load commuting_cost.mat;  
load calib_auxiliary.mat;

load fundamentals_1950.mat;
load calib_1950_auxiliary.mat;

load calib_qadjusted.mat; 
load noagglomeration.mat;         
load counterfactualpath.mat;     

Rdata=Rdata_all; Ldata=Ldata_all;  
totalpop=sum(Rdata(:,2:end),1);
N=size(ID,1);
NT=6; 
center=cbddist<2;
center=center*center';
% ==================== EXOGENOUS FUNDAMENTALS (50-75) =====================
a_mat=[a50, aa];   a_mat=a_mat./repmat(geomean(a_mat),N,1);
b_mat=[b50, bb];   b_mat=b_mat./repmat(geomean(b_mat),N,1);
A=a_mat.*(Ldata(:,2:end)./repmat(Area,1,NT)).^alphaest;  % overall productivity
B=b_mat.*(Rdata(:,2:end)./repmat(Area,1,NT)).^betaest;   % overall amenities
% ================= COMPUTE BASELINE COMMUTING PATTERN ====================
Lcom=zeros(N,N,NT+1); 
Lcomcenter=zeros(N,N,NT+1); 
Lcom(:,:,1)=Lcom45; 
Lcomcenter(:,:,1)=center.*Lcom45; 
Ymat=[Y50,Y55,Y60,Y65,Y70,Y75]; 
LQmat=zeros(N,NT);
LQmat(:,NT)=Qadj(:,NT);
for t=1:NT-1 
    tt=NT-t;
    theta=thetavec(tt); 
    if t < NT-1; LQmat(:,tt)=Qadj(:,tt).*(LQmat(:,tt+1).^(rho*(1-theta))); 
    else;        LQmat(:,tt)=Qadj(:,tt).*(LQmat(:,tt+1).^(rho50*(1-theta))); end
end 
for t=1:NT
    lambda_R=Kmat(:,:,t).*((repmat(LQmat(:,t),1,N).^(-gammaH/gammaL)).*(repmat(Ymat(:,t),1,N).^(1/gammaL))).^(rho/sigma);
    lambda_R=lambda_R./repmat(sum(lambda_R,1),N,1);
    if t==1; theta=theta50; 
    else     theta=thetavec(t-1); end 
    Rpost=Rdata(:,t+1); 
    Rpre=Rdata(:,t); 
    Lcom(:,:,t+1)=(1-theta).*Lcom(:,:,t)+lambda_R.*repmat((Rpost-(1-theta).*Rpre)',N,1);
    Lcomcenter(:,:,t+1)=center.*Lcom(:,:,t+1);
end
Lcomcenter_cf=repmat(center,1,1,NT+1).*Lcom_cf;

% =================== COMPUTE COMMUTING PATTERNS =============================
dLcom=zeros(N,N,NT);  
dLcom_cf=zeros(N,N,NT); 
for t=1:NT
    if t==1; theta=theta50; 
    else     theta=thetavec(t-1); end 
    dLcom(:,:,t)=Lcom(:,:,t+1)-(1-theta).*Lcom(:,:,t); 
    dLcom_cf(:,:,t)=Lcom_cf(:,:,t+1)-(1-theta).*Lcom_cf(:,:,t);
end

% ============== COMPUTE FLOW UTILITY =====================================
u_base=zeros(N,N,NT); 
u_cf=zeros(N,N,NT); 
for t=1:NT 
    w=(Qadj(:,t).^(-gammaH/gammaL)).*(A(:,t).^(1/gammaL));
    u_base(:,:,t)=repmat(B(:,t)',N,1).*repmat(w,1,N).*repmat(Qadj(:,t)',N,1).^(-mu)./Kappa_mat(:,:,t);
    w_cf=(Q_cf(:,t).^(-gammaH/gammaL)).*(A_cf(:,t).^(1/gammaL));
    u_cf(:,:,t)=repmat(B_cf(:,t)',N,1).*repmat(w_cf,1,N).*repmat(Q_cf(:,t)',N,1).^(-mu)./Kappa_mat(:,:,t);
end

% =============== COMPUTE WELFARE MEASURES ================================
% flow utility 
UU=zeros(N,N,NT); 
UU_cf=zeros(N,N,NT);  
UU(:,:,NT)=u_base(:,:,NT); 
UU_cf(:,:,NT)=u_cf(:,:,NT);  
for t=1:NT-1
    UU(:,:,NT-t)=u_base(:,:,NT-t).*(UU(:,:,NT-t+1).^rho);  
    UU_cf(:,:,NT-t)=u_cf(:,:,NT-t).*(UU_cf(:,:,NT-t+1).^rho); 
end

% commuting changes 
dL=zeros(N,N,NT); 
dL_cf=zeros(N,N,NT);  
dL(:,:,NT)=1; 
dL_cf(:,:,NT)=1; 
for t=1:NT-1
    dL(:,:,NT-t)=(dLcom(:,:,NT-t+1).^(-thetavec(NT-t))).*(dL(:,:,NT-t+1).^rho); 
    dL_cf(:,:,NT-t)=(dLcom_cf(:,:,NT-t+1).^(-thetavec(NT-t))).*(dL_cf(:,:,NT-t+1).^rho);
end
% last term (level term) 
M=zeros(NT,1); msize=1; M(NT)=1; 
for t=1:NT-1
    M(NT-t)= ((thetavec(NT-t).*msize).^thetavec(NT-t)).*(M(NT-t+1).^rho);
end

% welfare changes between counterfactuals and baseline (dV) 
dv_cf=zeros(NT,1);  
for t=1:NT
    ww=sum(sum((Lcom(:,:,t+1)./totalpop(t)).*(UU(:,:,t).^(rho/sigma)).*dL(:,:,t).^rho)); 
    ww_cf=sum(sum((Lcom_cf(:,:,t+1)./totalpop(t)).*(UU_cf(:,:,t).^(rho/sigma)).*dL_cf(:,:,t).^rho)); 
    dv_cf(t,1)=(sigma/rho).*log(ww_cf./ww);
end

% focus on center 
dvcenter_cf=zeros(NT,1); 
for t=1:NT
    ww=sum(sum((Lcomcenter(:,:,t+1)./sum(sum(Lcomcenter(:,:,t+1)))).*(UU(:,:,t).^(rho/sigma)).*dL(:,:,t).^rho)); 
    ww_cf=sum(sum((Lcomcenter_cf(:,:,t+1)./sum(sum(Lcomcenter_cf(:,:,t+1)))).*(UU_cf(:,:,t).^(rho/sigma)).*dL_cf(:,:,t).^rho)); 
    dvcenter_cf(t,1)=(sigma/rho).*log(ww_cf./ww);  
end

% welfare dynamics (V)
v_base=zeros(NT,2); 
v_cf=zeros(NT,2); 
for t=1:NT
    v_base(t,1)=sigma.*log(sum(sum((Lcom(:,:,t+1)./totalpop(t)).*(UU(:,:,t).^(rho/sigma)).*dL(:,:,t).^rho)))+sigma.*log(M(t).^rho);  
    v_base(t,2)=sigma.*log(sum(sum((Lcomcenter(:,:,t+1)./sum(sum(Lcomcenter(:,:,t+1)))).*(UU(:,:,t).^(rho/sigma)).*dL(:,:,t).^rho)))+sigma.*log(M(t).^rho);
    v_cf(t,1)=sigma.*log(sum(sum((Lcom_cf(:,:,t+1)./totalpop(t)).*(UU_cf(:,:,t).^(rho/sigma)).*dL_cf(:,:,t).^rho)))+sigma.*log(M(t).^rho);  
    v_cf(t,2)=sigma.*log(sum(sum((Lcomcenter_cf(:,:,t+1)./sum(sum(Lcomcenter_cf(:,:,t+1)))).*(UU_cf(:,:,t).^(rho/sigma)).*dL_cf(:,:,t).^rho)))+sigma.*log(M(t).^rho);
end


result=array2table([round(dv_cf,3),round(dvcenter_cf,3)],'VariableNames',{'welfare_cf',...
    'centerwelfare_cf'}); 
writetable(result,'../../output/quant/output_welfare_counterfactualpath.csv');

